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Abstract 

We propose a computationally efficient technique for extrapolating 
seismic waves in an arbitrary isotropic elastic medium. The method 
is based on factorizing the full elastic wave equation into a product of 
pseudo-differential operators. The method extrapolates displacement fields, 
hence can be used for modeling both pressure and shear waves. The pro- 
posed method can achieve a significant reduction in the cost of elastic 
modeling compared to the currently prevalent time- and frequency-domain 
numeric modeling methods and can contribute to making multicomponent 
elastic modeling part of the standard seismic processing work flow. 

1 Introduction 

Extrapolation of seismic wave fields in depth using one-way propagation opera- 
tors is an efficient alternative to time- and frequency-domain modeling with the 
full wave equation, particularly in seismic migration applications (see [2], PQ). 
While one-way extrapolators have long been established as key components of 
the seismic imaging toolbox for isotropic acoustic media, extrapolation of elastic 
wave fields is still carried out by solving the full elastodynamic system either in 
the time or frequency domain, either approach being computationally expensive. 
The high computational cost of wave extrapolation in elastic media is one of the 
barriers to a widespread adoption of multicomponent seismic in industrial appli- 
cations. Some progress has been made recently in the development of efficient 
one-way methods for certain simplest anisotropic elastic models (e.g., vertically 
transversally isotropic or tilted transversally isotropic media - see [TO], [S], [E]) 
However, these methods use the "pseudoacoustic" approximation (see [3]) and 
are used for a kinematically accurate propagation of pressure waves only. 

In this paper we present a method for one-way frequency-domain extrapo- 
lation of displacement fields in an elastic isotropic medium. The approach of 
this paper is based on factorizing clastic wave equation using pseudo-differential 
operators without introducing stress-related unknown functions into the equa- 
tions. Our approach is conceptually similar to the derivation of the acoustic 



single square-root equation (see 0) except the resulting factorized propagation 
operators can not be obtained analytically but are computed numerically. 



2 The Method 

We start with the wave equation governing the displacement in an arbitrary 
heterogenous isotropic elastic medium in the Navier form (see [9]): 
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where u l denote the components of a displacement field, p is the shear mod- 
ulus, v is Poisson's ratio for the medium, and p is the density. In this paper 
we consider a heterogenous elastic medium under the assumption of local ho- 
mogeneity - otherwise the elastic moduli would not be factored outside of the 
differentiation operators in equation [l] However, our method can be extended 
to cover the case when local homogeneity assumption is dropped. "Freezing" 
the coefficients of equation [T] and applying the Fourier transform in time and 
horizontal variables x\ — x,X2 — y, and substituting 
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where A is the Lame coefficient (see [Z],[9]), we get 
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where k x , k y are horizontal wave numbers and u) is the frequency. The left-hand 
side of system [3] is the result of an ordinary differential operator applied to a 
vector-function u = (it 1 ,u 2 ,w 3 ) and parametrized by horizontal wave numbers. 
In the present form equations [3] cannot be used for computationally efficient ex- 
plicit depth extrapolation in a heterogeneous medium; however, these equations 
can be used for modeling displacements by solving a series of boundary-value 
problems (see [S]). In [S] it was suggested that equations [3] might be factorized 
in such a way as to allow solving them by alternating one-way extrapolation 
in opposite directions. More specifically, we seek a factorization of operator 
equation [3] of the form: 
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and A, B are 3x3 matrices with components that are complex- valued functions 
of the horizontal wave numbers, I is the 3x3 identity matrix. Performing the 
multiplication in equation [4] and using equation [3j we obtain: 



A(k Xl k y )B(k Xl k y ) + c CJ [j4(fc a; , k y ) + B(k x , k y )] — P(k x , k y ), 
A(k x ,k y )E(X, n) + E(X, n)B(k x ,k y ) + 2c M E(X, p) = S{k x , k y ), 
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Combining equations [6] and [7J we get the following equation for the operators 
A and B : 

A(k x , ky)B(k x , k y ) + c LJ (A(k x , k y ) + B(k xi k y )) = P(k x , k y ), 

E(X, n)B{k x , k y ) + A{k x , k y )E(X, n) = S{k x , ky), (8) 



where 



S(k x ,k y ) = S(k x ,k y ) - 2c u E(X,fj,). 



(9) 



Equations |4j [8] in combination with equations [7] and [9] suggest the following 
procedure for extrapolating solutions to system [T] in depth: 

1. Solve the system of matrix equations [8] for A, B, for each pair of horizontal 
wave numbers k x , k y and two reference values of each elastic parameter 
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2. Evaluate 
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from the initial conditions and assign the value to an auxiliary function 
u(x,y,z = 0); 
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3. Solve 

(e(X, + A(-id x , -id y ) + cj\ u(x, y,z)=0 (10) 

by downward continuing in depth, using the formula 

u(x, y,z + Az) — exp [— AzE (A(— id x , —id y ) + c^T)] u(x, y, z). (11) 

4. Perform each step of the depth extrapolation for four combinations of the 
reference elastic parameters, then apply inverse Fourier transform to the 
four fields and interpolate at each spatial point of the depth slice using 
true X(x,y), /i(x,y) as e.g. in the PSPI method (see pQ). 

5. After reaching the desired maximum depth, find the solution u by upward 
extrapolation: 

d_ 

' dz 



E(\,n)— + B(-id X) -idy) + c u I ) u(x,y,z) = u(x,y,z). (12) 



6. Repeat the above steps for each frequency component u(uj,x,y,z). 
The above algorithm is stable if the spectrum of matrix 

A{k x ,ky)+cJ (13) 
is not in the interior of the left half-plane, and the spectrum of 

B(k x ,k y )+c u I (14) 

is not in the interior of the right half-plane. While the above algorithm tries to 
mimic the two-way wave propagation, it is effectively just an approximation to 
the propagation process as it ignores the interaction between the up and down- 
going wave at intermediate depth steps. A less accurate alternative would be 



to downward-continue the wave field using equation 11 in a way similar to the 
one-way depth extrapolation using the scalar square- root equation (see [2J , PQ ) - 
The latter approach would be unable to image any dips beyond 90°, however, 
it would reduce the cost of extrapolation by a further factor of 2. Note the cost 
of solving equation [TU] in depth is roughly three times that of solving the scalar 
square-root equation. 

The above analysis may be extended to the case of an arbitrary anisotropic 
elastic medium. The fact that the components of the pseudo-differential opera- 
tor matrices A(—id x ,—id y ),B(—id x ,—id y ) are not given in an analytical form 
but are only computed numerically does not limit their applicability. 

Factorization of system [3 in the elastostatic case was one of the approaches 
mentioned by the author in 5J. However, the one-way extrapolation technique 
is mostly useful for elastodynamic problems as the passband of the factorized 



depth extrapolators (e.g., as in equation 11) narrows down to zero with the 
temporal frequency passing to the zero static limit. Note that equation [T] uses 
elastic parameterization that degenerates into a singularity if the shear modulus 
is equal to zero. This is not causing any problems with purely acoustic wave 
extrapolation as the singularity is effectively removed from equations [3] by the 
substitution 
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Figure 1: The phase of a phase-shift operator corresponding to the maximum 
imaginary part of the eigenvalues of operator 15 Multicomponent "phase-shift" 



is defined by three such scalar phase-shift operators and a 3 x 3 matrix Q of 
equation |16| 



3 Implementation and Results 

The system of matrix equations [8] is solved only once for each triple of temporal 
frequency and elastic moduli values, and for each pair of horizontal wave num- 
bers. In our prototype implementation of the one-way extrapolator we compute 
the matrices A, B at the beginning of the frequency loop and subsequently use 
the tabulated matrices in the depth extrapolation loop (inside the frequency 
loop). A more efficient approach can be employed in a production implemen- 
tation of the extrapolation method: system [8] can be solved using Newton's 
method (see e.g. [1]) in a one-off computation for each set of the temporal 
frequency, elastic moduli and horizontal wave numbers and stored in a look- 



up table. The symmetry of the extrapolation operators 11 that appear to be 



multi-component counterparts of the acoustic phase-shift operator (see [2]) can 



5 



be exploited to achieve a substantial reduction in the size of the precomputed 
operator tables. Figure [T] is the plot of the maximum of the imaginary parts of 
the three eigenvalues of operator 



K 
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(15) 



within its passband. 
of figures |2|3|4|5 



The operator is the one used later to produce images 

are zero 



The real parts of the eigenvalues of operator 15 



within the operator passband and negative outside. The imaginary parts of 
the other two eigenvalues exhibit similar behavior. Operator K of equation 
[TB] is the logarithm of the extrapolation operator [TT] and the spectral plot of 
figure [T] corresponds to the phase of the phase-shift extrapolator in the acoustic 
case (see [I]). The crucial difference in the elastic multicomponent case is that 
the multicomponent "phase-shift" is defined by three such scalar phase-shift 
operators with phases 4>\ , <p2 , 4>3 , and a unitary operator Q, determined by the 
eigenvector expansion of K as follows: 



K = Q 



i<t>i 
i4>2 
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(16) 



The pass bands of the three phase shift operators are, generally, different, 



but the real parts of the eigenvalues of 15 are non-positive across all three pass 
bands. 

Figures |2|3|4|5| demonstrate the result of applying our method to extrapo- 
lating displacement waves from a concentrated impulse at the surface. Medium 
parameters used in this test were 316 m/ s shear- wave velocity 

us = VmTp 

and 632 m/ s pressure- wave velocity 

Up = y/(X + 2(l)/ P . 

The extrapolation grid was 128 x 128 x 128 with a 5 m step, frequency range 
2-12 Hz with 1 Hz step. The values of the elastic parameters used in this test 
are uncharacteristically low for seismic applications and were chosen solely for 
the purpose of fast small-scale simulation on a single-core PC using Matlab. 

Since the impulse at the surface is an asymmetric horizontal displacement 
but can be assumed to be symmetric in the vertical direction, our waves are ef- 
fectively a combination pressure and shear waves for the horizontal components, 
while the vertical displacement wave should kinematically match the pressure 
wave. And indeed, the pressure wave plot [2] and vertical displacement plot [3] 
exhibit excellent kinematic agreement. 

The horizontal wave component plots [4] and [5j on the other hand, show 
pressure- and shear-wave images, both correctly positioned in agreement with 
the velocity values used in the simulation. Boundary reflections and low fre- 
quency content cause some imaging artifacts that are not related to the method. 
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Figure 2: Pressure wave extrapolated from an impulse displacement, 2-12 Hz, 
128 x 128 x 128 grid, 5 m step, 316 m/s shear-wave and 632 m/s pressure-wave 
velocity. 

4 Conclusions and Discussion 

The method presented in this paper can be used in seismic migration algorithms 
in order to achieve a substantial reduction of run time in comparison with the 
reverse time migration. More specifically, stability of the time-domain modeling 
typically utilized in the reverse-time migration requires time steps significantly 
smaller than the time resolution of seismic data (see (T]). Depth extarpola- 
tion of wave fields using one-way equations [10] and [12] can be performed for an 
arbitrary frequency range. Extrapolating wave fields in the frequency domain 
using two-way system [3] would require solving a large sparse system of equations 
using e.g. finite elements method, still posing significant computational chal- 
lenges for inhomogeneous media. However, the one-way extrapolation method, 
while limited in dip and less accurate in terms of amplitudes, lends itself to 
efficient implementation using e.g. PSPI or finite differencing. Furthermore, 
the approach of this paper can be expected to apply to more complex elastic 
anisotropic models (see [3]) and may be developed into a computationally effi- 
cient alternative to the existing pseudo-acoustic anisotropic modeling methods 
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Figure 3: Vertical component of a wave extrapolated from an impulse displace- 
ment, inline section, 2-12 Hz, 128 x 128 grid, 5 m step, 316 m/s shear-wave and 
632 m/s pressure- wave velocity. 

while allowing easy separation of pressure and shear waves. 
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